Spectra of random graphs with arbitrary expected degrees 
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We study random graphs with arbitrary distributions of expected degree and derive expressions 
for the spectra of their adjacency and modularity matrices. We give a complete prescription for 
calculating the spectra that is exact in the limit of large network size and large vertex degrees. We 
also study the effect on the spectra of hubs in the network, vertices of unusually high degree, and 
show that these produce isolated eigenvalues outside the main spectral band, akin to impurity states 
in condensed matter systems, with accompanying eigenvectors that are strongly localized around 
the hubs. We also give numerical results that confirm our analytic expressions. 
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I. INTRODUCTION 

The topology of complex networks, such as social, bi- 
ological, and technological networks, can be represented 
in matrix form using an adjacency matrix or any of sev- 
eral other related matrices such as the graph Laplacian 
or the modularity matrix P, Q ■ The spectral properties 
of these matrices — their eigenvalues and eigenvectors — 
are related to a range of network features of scientific 
interest, including optimal partitions [3, 0] , percolation 
properties Q, community structure an d the behav- 

ior of network dynamical processes such as random walks, 
current flow, diffusion, and synchronization @, Q . As a 
result, the study of network spectra has been the subject 
of considerable research effort for some years. This effort 
has taken a number of forms. One has been the study 
of the spectra of empirically observed networks, which 
can be calculated by numerical means for networks of 
size up to hundreds of thousands of vertices [8j, [lfj . An- 
other, which is the topic of this paper, is the study of the 
spectra of model networks. A fundamental question we 
would like to answer is how particular structural features 
of networks are reflected in network spectra, and model 
networks provide an ideal setting in which to investigate 
this question. 

Some results in this area have been known for a long 
time. For example, the very simplest of network mod- 
els, the Poisson random graph, studied as far back as 
the 1950s by Erdos, Renyi, and others [III G3: nas a 
symmetric adjacency matrix whose elements are inde- 
pendent identically-distributed random variables. Such 
matrices are known, subject to some conditions but re- 
gardless of the precise distribution of their elements, to 
have a universal spectrum obeying the Wigner semicircle 
law, and eigenvectors that are distributed isotropically at 
random, except for the leading eigenvalue and eigenvec- 
tor, whose values are governed by the Perron-Frobenius 
theorem and the average degree of the network [l3l-|20|. 

As we have come to understand in the last decade, 
however, the random graph is a poor model for the struc- 
ture of real- world networks. In particular, the frequency 
distribution of the degrees of vertices in the random 
graph is Poissonian, while the degree distribution of most 
real-world networks is highly right-skewed, often having 



a power-law or exponential tail of "hubs" with degree 
far above the mean. Luckily, it turns out to be possi- 
ble to create generalizations of the basic random graph 
that incorporate arbitrary degree distributions, including 
skewed distributions, the best-known such model being 
the so-called configuration model [2l|, [22[. The config- 
uration model is solvable exactly for many of its struc- 
tural properties, including its complete component struc- 
ture (22| - |23 | and percolation properties [2^, [26[ , and the 
results have led us to a better understanding of the pro- 
found effect the degree distribution has on network topol- 
ogy- 

In this paper we study the spectral properties of the 
configuration model. Motivated by recent developments 
in random matrix theory, we derive a simple recipe for 
calculating the spectrum of the adjacency matrix of the 
model. We show that the spectrum is composed of three 
fundamental elements, all of which have clear correlates 
in the structure of the network. The elements are: (1) the 
leading eigenvalue, which is dictated primarily by the av- 
erage network degree; (2) a continuous band or "bulk 
spectrum," analogous to the Wigner semicircle but tak- 
ing a different shape; and (3) in some but not all cases, 
additional eigenvalues outside of the continuous band 
which correspond to the hubs in the network and which 
have eigenvectors that are strongly localized about those 
hubs. 

In addition to our analytic developments, we also con- 
firm the form and behavior of each of these elements with 
numerical calculations on example networks generated 
using the configuration model. 

A number of previous authors have examined 
the spectral properties of the configuration model. 
Farkas et al. |10| performed numerical calculations on 
large samples generated using the model and demon- 
strated that there are clear deviations from the semicircle 
law for non-Poisson choices of the degree distribution, 
and especially for power-law distributions. Dorogov- 
tsev et al. [27| gave an analytic route to the full spectrum, 
though their method is complex, involving the solution 
of a nonlinear integral equation containing Bcssel func- 
tions, which at present can only be done approximately. 
Chung et al. |28[ gave a rigorous derivation of the ex- 
pected value of the largest eigenvalue in the spectrum in 
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the limit of a dense network. Our calculations extend 
these studies by providing a simple derivation of the full 
spectrum which is exact in the limit of large vertex de- 
grees and confirms earlier findings while shedding new 
light on features of the spectrum and their implications 
for network structure. 



II. THE MODEL 

In this paper we study the spectral properties of the 
configuration model — or, more precisely, a slight variant 
of the model, as we now describe. 

The configuration model is a model of an undirected 
random graph with a specified number of vertices n and a 
given degree sequence. In this model one first specifies a 
degree sequence, meaning one specifies the degree of each 
of the n vertices. Let the degree of vertex i be denoted ki 
and let us visualize the degree as ki ends or "stubs" of 
edges emerging from the vertex. Then the configuration 
model is defined as the ensemble of pairwise matchings of 
stubs in which every matching appears with equal prob- 
ability. That is, a configuration model network with the 
given degree sequence is generated by repeatedly choos- 
ing two stubs uniformly at random from those available 
and joining them together to form a complete edge. This 
process continues until all stubs have been joined and no 
unattached stubs remain. (For this to work, the number 
of stubs must be even, and hence the model is defined 
only for degree sequences whose sum J^. ki is even.) 

The configuration model provides a way to generate 
networks that have any degree sequence we desire while 
being essentially random in other respects — there are no 
correlations or long-range structure in the configuration 
model ensemble. 

A crucial feature of the configuration model for our 
purposes will be the expected number of edges between 
a vertex pair. It is straightforward to show, given the 
degree sequence, that the expected number of edges be- 
tween vertices i and j is equal to kikj/2m in the limit 
of large network size, where m = ^ ki is the num- 
ber of edges in the network. Note that it is possible 
to generate networks with multi-edges — pairs of vertices 
connected by more than one parallel edge. The actual 
number of edges between vertices i and j is multinomi- 
ally distributed with mean kikjjlm. 

However, edges in the configuration model are not sta- 
tistically independent. Since the degrees of vertices are 
fixed, the presence of an edge from vertex i to vertex j 
makes it less likely that there will be an edge from i 
to any other vertex, and hence edges that share a com- 
mon end are correlated. When degree is large the cor- 
relations become small and the multinomial distribution 
of edge number becomes approximately Poisson, but for 
networks with finite average degree the correlations will 
always be present and may be significant. 

These correlations make analysis of the model more dif- 
ficult and so in this paper we consider a modified model 



in which the number of edges between each pair of ver- 
tices is defined to be an independent random variable 
with mean kikj/2m and value drawn from a Poisson dis- 
tribution with that mean. In this model, ki becomes the 
expected degree of vertex i and m is the expected total 
number of edges. When degrees become large, which is 
the primary regime that we consider in this paper, the 
actual degrees will be narrowly peaked about their ex- 
pected values, so the properties of the variant model and 
the standard configuration model, including the spectral 
properties that we study, become the same. This model 
(or slight variants of it) has been studied previously by 
a number of authors, notably Chung and Lu jlgj], with 
whose work it is perhaps most strongly associated. 

In this paper we consider networks in the limit of large 
size n with expected vertex degrees drawn from a fixed 
probability density p(k), so that p{k) dk is the fraction 
of vertices with expected degree in the interval from k 
to k + dk. (Note that expected degree need not be an 
integer, although one is free to choose it to have integer 
values if one wishes.) More precisely, we consider a se- 
quence of networks of increasing size with fixed expected 
degrees ki and additional degrees drawn from p(k) as n 
becomes larger. Thus for finite n the expected degree of 
any particular vertex i remains constant as n becomes 
large and the empirical degree distribution converges to 
p(k) in the large-n limit. 

The adjacency matrix A of a network generated ac- 
cording to this model is the nxn symmetric matrix with 
integer elements Ay equal to the number of edges be- 
tween vertices i and j. Our primary goal in this paper is 
to calculate the average spectrum of the adjacency matrix 
within the model ensemble, which we do in two stages. 
We write the matrix as 

A=(A)+B, (1) 

where (A) is the ensemble average of A, which has ele- 
ments (Aij) = kikj/2m, and B is the deviation from that 
average. Our approach is first to calculate the spectrum 
of the matrix B, whose elements are, by definition, in- 
dependent random variables with zero mean, although 
crucially they are are not identically distributed. Once 
we have the spectrum of B then the spectrum of A is 
calculated from it in a separate step. 

The matrix B is of interest in its own right. It has 
elements 

Bij = - (Aij) = Aij (2) 

This matrix is known as the modularity matrix, and 
forms the basis for one of the most widely used methods 
for detecting modules or communities in networks 0, 0] ■ 
The methods described in this paper thus give us the 
spectra of both the adjacency matrix and the modularity 
matrix. 

Note that the elements of the modularity matrix have 
variance the same as the elements of the adjacency matrix 
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which, since they are Poisson distributed, have variance 
equal to their mean fc i fc J /2m. Hence 



which will be important shortly. 



(3) 



III. SPECTRUM OF THE MODULARITY 
MATRIX 

As discussed in the previous section, we will first cal- 
culate the spectrum of the modularity matrix B, defined 
by Eq. ([2]), then calculate the spectrum of the adjacency 
matrix from it in a separate step. We begin by developing 
some fundamental notions concerning random variables 
that will be important for our derivations. 

Suppose we have two independent random variables, x 
and y 7 ordinary scalar variables, with probability densi- 
ties p x (x) and p y (y)- What is the probability that their 
sum x + y will have a particular value z? The answer to 
this question is well known and simple. The probability 
density for z is 



p(z) 



p x {x)p y (y)5(x + y - z) dxdy 



= / p x {x)p y {z - x) Ax 



(4) 



which is the convolution of the two distributions. Simi- 
larly we can ask for the probability that the product xy 
has value z, which is given by the multiplicative convo- 
lution 



p{z) 



Px(x)Py(y)5(xy - z) dxdy 

P x (x)p y (z/X) — . 



(5) 



A scalar random variable can be thought of as the sin- 
gle eigenvalue of a 1 x 1 random matrix. A 1 x 1 matrix is 
diagonal by definition and its one eigenvalue is trivially 
equal to its one element. A natural generalization of the 
convolution results above is to ask what their equivalent 
is for larger random matrices, 2 x 2, 3 x 3, and so forth, 
where we will confine ourselves to symmetric matrices, 
so that the eigenvalues are real. That is, if we know 
the probability density of the eigenvalues — the so-called 
spectral density — of two independent symmetric random 
matrices, what is the spectral density of their sum or 
product? The answer is no longer a simple convolution, 
because matrices do not in general commute, so what is 
the appropriate generalization? Unfortunately, this ques- 
tion does not have a straightforward answer because it 
turns out that a knowledge of the spectral densities alone 
is not enough. In general one needs to know the distri- 
bution of the entire matrices to calculate the spectral 
density of their sum or product. There is, however, one 
case in which relatively simple results apply, which is 



when the eigenvectors of the two matrices are themselves 
random and uncorrelated. 

Recall that the eigenvectors of a symmetric matrix 
are orthogonal — for an n x n matrix they define a set 
of orthogonal axes in an rt-dimensional vector space. 
Thus if we have two random symmetric matrices, the 
eigenvectors of one can always be transformed into the 
eigenvectors of the other by a suitable rotation and/or 
reflection — in other words by a suitable unitary trans- 
formation. If for different choices of the random ma- 
trices the transformations needed to do this are dis- 
tributed isotropically — if all possible such transforma- 
tions are equally likely — then the random matrices are 
said to be free. Loosely, one can say that two random 
matrices are free if the angle between their eigenvectors 
is also random. The mathematics of free random vari- 
ables has been developed extensively since the 1990s and 
is known by the name of free probability theory [30j . 

The crucial observation now is the following: for free 
matrices the spectral density of their sum or product is a 
function only of the individual spectral densities. It turns 
out that one no longer needs to know the entire distribu- 
tion of the matrices themselves and well-defined general- 
izations of the convolution equations, Eqs. (|4|) and (|5|), 
exist. For the sum of two matrices the appropriate gener- 
alization is known as the free convolution or free additive 
convolution] for the product of matrices it is the free 
multiplicative convolution. Thus if two symmetric ran- 
dom matrices have spectral densities p x {x) and p y (y), 
then the spectral density of their product is the free mul- 
tiplicative convolution 



P{Z) = (p X Mp y )(z), 



(6) 



where Kl denotes the convolution. Although this defines 
the convolution in principle, it does not tell us how to 
calculate it. We will come to that in a moment, but first 
let us return to the configuration model and see why this 
is a useful result. 

We wish to calculate the spectral density of the mod- 
ularity matrix B, which for an undirected network is 
a symmetric random matrix whose elements have zero 
mean but different variances, equal to kikj/2m — see 
Eq. ([3]) . Let us define a normalized modularity matrix B 
by 



B = D 1/2 BD 1/2 



(7) 



where D is the diagonal matrix with elements ki . B has 
elements Bij = B^ / >J ki kj , so that each is divided by a 
factor proportional to its standard deviation and hence, 
though not identically distributed, all elements now have 
the same variance, equal to l/2m. So long as the vertex 
degrees are large, matrices with this property are known 
to have an eigenvector basis oriented isotropically at ran- 
dom and to have spectral density obeying the Wigncr 
semicircle law [13j-|20(, which for our particular matrix 
takes the form 



o c (z) = ^-^lc~ 



c 2 z 2 , 



(8) 
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where c = 2m I n is the average degree in the network. 
The requirement that vertex degrees be large is necessary 
because deviations from the semicircle law are known to 
arise for very sparse matrices [3l|. For small degrees, 
therefore, the results given here will only be approximate. 

Now consider an eigenvalue z of the modularity ma- 
trix B itself, satisfying Bb = zh where b is the cor- 
responding eigenvector. Multiplying by D 1 / 2 , writing 
B = D^BD 1 / 2 , and defining v = D 1 / 2 b, this can also 
be written 



DBv 



(9) 



In other words the modularity matrix has the same eigen- 
values as the matrix DB, which is the product of the di- 
agonal matrix D, which by definition has spectral density 
equal to the degree distribution p(k), and the symmetric 
matrix B, with spectral density p c {z) given by Eq. jSJ. 

But it is precisely to the products of such random ma- 
trices that Eq. relates, and hence, applying that equa- 
tion, we arrive at the principal result of this paper: the 
spectral density of the modularity matrix for a network 
with arbitrary expected degrees is equal to the free mul- 
tiplicative convolution of the degree distribution with the 
Wigner semicircle. That is, the spectral density p(z) is 
given by 



p(z) = (pBp c )(z), 



(10) 



where p{k) is the distribution of expected degrees and 
p c (z) is given by Eq. (HJ). 

This result is of immediate practical utility. Numerical 
methods exist for c omp uting free multiplicative convolu- 
tions efficiently [32 , |33j , which means we can use existing 
numerical packages to compute spectral densities easily 
and rapidly for a wide range of degree distributions. 

For the purposes of the present paper, however, we 
would like to know more. In particular, we would like 
explicit formulas for calculating the spectral density in 
the general case. Unfortunately, the free multiplicative 
convolution has no simple expression for matrices of fi- 
nite size, but in the limit of large size — which is also the 
limit of a large network — suitable expressions do exist. 
Specifically, for a spectral density p we can define a func- 
tion 



x p(x) dx 



(11) 



which is called the Cauchy transform of xp(x). Then 
if p is the free multiplicative convolution of two other 
distributions p and p c as in Eq. (|10[) . it can be shown 
that 

r^(u) = ^-T- 1 (u)T^(u), (12) 
' u + 1 F ' 

where T _1 denotes the functional inverse of T, and T p 
and r Pc are defined by analogy with (|TT|) : 



r P (*) 



x p c (x) dx 



(13) 



Substituting Eq. ((SJl into the second of these, we have 



Z7T 



:V4c 



dx 



2/>/c Z — X 

±cz(z±y/z 2 -4/c) 



1. 



(14) 



The ambiguity in the sign of the square root arises be- 
cause of a branch cut in the evaluation of the integral, 
but it can be shown that the final result for the free con- 



volution never depends on the choice of sign [34| . Here 
we take the negative sign, since it makes some of the fol- 
lowing steps cleaner. Rearranging for z as a function of 
r„ we then find that the functional inverse is 



(15) 

and substituting into (fT2|) we get 

T- 1 (u) = ^T-\u). (16) 

Evaluating this equation at the point u = T p (z) gives 
z = y/T p (z)/c T p 1 (T p (z)), which can be rearranged to 
read 



T p (z) = T p (z^c/T p (;. 



(17) 



For convenience we define h(z) = ^jT p {z)/c and Eq. (fl~7|) 
becomes 



ch 2 (z) =T p (z/h(z)) = 
or, more simply, 



kp{k) dk 



h(z) = - 

c Jo 



z/h(z) — k ' 

kp(k) dk 
z — kh(z) 



(18) 



(19) 



If we can solve this equation for h(z) then the Cauchy 
transform of xp(x), Eq. (|11[) , is given by T p (z) = ch 2 (z). 
To recover p itself from the Cauchy transform we note 
that for real x and rj 



— Im — 

7T x 



1 



17] 



.2 _|_ ^2 



(20) 



which is a Lorentzian of width ry and area 1, and hence 
in the limit as rj — > + becomes equal to a delta-function: 



— lim Im 

7T 7/->-0+ X + 



177 



5(x). 



(21) 



Thus 



zp(z) = / xp(x)S(z — x) dx 



lim Im 

7T ^->-0+ 



xp{x) 



dx 



177 



= lim Imr„(z + ir?). 

7T t?^0+ ' ' 



(22) 
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This is the Stieltjes-Perron inversion formula. Setting 
T p (z) = ch 2 (z) it tells us that the spectral density of the 
configuration model is given by 



p{z) 



Im h 2 {z), 



(23) 



where the imaginary part is taken in the limit as z tends 
to the real line from above. 

Equations and (|2"3"|) give us a complete recipe for 
calculating the spectrum of the modularity matrix. We 
note that equations equivalent to these have been derived 
in other contexts in the literature on random matrices. 
See for example the results on band matrices in Refs. [3E 

m. 
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A. Example solutions 



The solution of Eqs. (fTi?]) and ([2"B"f relies on our being 
able to compute the integral in Eq. (|T9| , whose difficulty 
depends on the particular choice of degree distribution. 
To give an example where the calculation is straightfor- 
ward, consider the standard Poisson random graph, for 
which all vertices have the same expected degree c and 
hence p{k) = S(k — c). Substituting into (TJ5]) and solving 
the resulting quadratic equation gives 



h(z) 



4c 



2c 



so that the spectral density is 



V4c 



2irc 



(24) 



(25) 



which recovers the standard semicircle distribution for 
the random graph. 

As a more general example, consider any distribution 
where the degrees take a set of £ discrete values d r , as 
they do for any integer- valued degree distribution of the 
type commonly considered for network models. Then 
p(k) = ^2 r —iPr^(k — d r ), where the coefficients p r sat- 
isfy J2 r Pr = !■ Then, from Eq. (fT9|) . 



h(z) 



J2 r =i Prd r /[z - d r h(z)] 

Sr=l Prd r 



(26) 



where we have used c = ^2 r p r d r - Thus h(z) is the root 
of a polynomial of degree £+1. For instance, if there are 
two discrete values of the expected degree, then 



h(z) = 



1 



Pidi + p 2 d 2 



Pidi 



Pidi 



z — d\h{z) z — d 2 h(z) 
which can be rearranged to give the cubic equation 

d\d 2 



(27) 



did 2 h 3 -(d 1 +d 2 )zh 2 



Pidi +p 2 d; 



■+z' 



h-z = 0. (28) 



FIG. 1: The spectral density p(z) for the degree distribution 
described in the text, in which vertices have expected degree 
d\ = 50 with probability pi = j and di = 100 with proba- 



bility P2 



The curve gives the analytic solution, derived 



from Eqs. (ffifl) and the histogram shows the results of 

numerical calculations for actual networks with n = 10 000 
vertices, averaged over 100 different networks. 



Of the three solutions to this equation one is always real, 
and hence (in light of Eq. (|2"3")l ) cannot give the spectral 
density. The remaining two are complex conjugates and 
so give results that differ only in sign, the positive sign 
being the one we are looking for. 

Figure [T] shows a plot of the resulting spectral density, 
Eq. (f2"3"|) . as a function of z for the case d\ = 50, d 2 = 100, 
and pi = 1 — p 2 = -7. The figure shows strong depar- 
ture from the semicircle law. Also shown are the results 
of direct numerical calculations of the spectra of simu- 
lated networks with the same degree distribution and the 
agreement between the analytic and numerical results is 
good. 



B. Features of the spectral density 

We can invoke additional properties of the free con- 
volution to better understand the spectrum of the mod- 
ularity matrix. Consider, for instance, the case where 
the expected degree distribution p{k) has compact sup- 
port, meaning that there are hard upper and lower limits 
to the expected degree a vertex may have. (The lower 
limit is trivial, since degrees must be non-negative, but 
the upper limit is not.) Since the semicircle distribution, 
Eq. ([8]), also has compact support, the spectral density 
of the modularity matrix is then a convolution of two 
compact distributions. In this scenario it can be shown 
that the bulk spectrum of the modularity matrix will also 
have compact support (3^|. Furthermore, given this ob- 
servation we can show that the spectrum will generically 
exhibit a sharp square-root decay at its edges. To see 
this, note that the central function h(z) in our theory is 
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the solution for h of an equation of the form f(h, z) = 
where z is given and 

\ 1 f°° kp(k)dk , 



(See Eq. (HI]).) From Eq. (J23J) we know that h is com- 
plex within the spectral band and real outside it and 
hence the edge of the band is the point at which complex 
solutions to f(h, z) = disappear. For analytic f(h, z) 
such a disappearance corresponds to the point at which 
an extremum of / with respect to h crosses the zero line. 
Denoting this point by (h, z) = (a, b) and performing an 
expansion about it to leading order in both h and z, we 
then have 



f^ z)= d J- {z -b) + ^L{h-af + 



(30) 



the terms in f(a,b) and df/dh vanishing at the ex- 
tremum. In the limit as we approach the band edge, 
therefore, the equation f(h, z) = takes the form 



g-z {z -V + W (h - a) 



0, 



(31) 



and hence, within the band, we have h(z) = a + iB^/b — z 
for some real constant B. Then the spectral density, 
Eq. (J221), is 



p{z) = C 



(32) 



where C is another real constant. A similar argument 
implies square-root behavior at the lower edge of the 
spectrum as well. The square-root form can be seen, for 
example, in the vertical sides of the spectrum in Fig. [1] 
We can also calculate the behavior of h(z) as z — >• b 
from above, for which Eq. (J3TJ) implies 



h(z) 



BVz - b, 



(33) 



with the same real constant B as before. Note that this 
implies that the limiting value of h(z) at the band edge 
is generically finite, but that the slope dh/dz diverges. 
This has important consequences for "hub" vertices — 
those with unusually high degree — whose effect on the 
spectrum displays a phase transition behavior that de- 
pends crucially on the functional form of h(z). We dis- 
cuss hub vertices in detail in Section fVTl 

These results apply for the case where the expected 
degree distribution is bounded. In cases where it is not 
we expected the spectral density of the modularity ma- 
trix to be similarly unbounded, having no band edge 
and generically inheriting the worst-case tail behavior 
of p(k). Similar observations have been made previ- 
ously by Chung et al. 28j] for a different matrix, the 
graph Laplacian. They note that a normalized version 
of the Laplacian, akin to our normalized modularity ma- 
trix B, should display a semicircle distribution, but that 
the Laplacian itself should have a spectrum that inherits 
the tail behavior of the degree distribution. 



IV. THE RESOLVENT AND THE STIELTJES 
TRANSFORM 

In the previous section we calculated the spectral den- 
sity of the modularity matrix for the configuration model. 
It is possible to calculate many other properties of the 
spectrum as well, as we now show. Our starting point 
for these calculations is the so-called resolvent matrix, 
which is the matrix function R(z) = (zl— B) _1 , where I 
is the identity. As we will see, a knowledge of the ensem- 
ble average of the resolvent allows us to calculate many 
things, including the spectral density of the adjacency 
matrix, the leading eigenvalue of the adjacency matrix, 
and the effect on the spectrum of network hubs. 

It also gives us an alternative, though perhaps less ele- 
gant, derivation of the results for the modularity matrix 
in the previous section. The spectral density p(z) of the 
modularity matrix can be defined as 



P{z) 



1 ™ 

71 ^ 



A-s), 



(34) 



where Xi are the eigenvalues of the matrix. Substituting 
for the delta- function from Eq. (|2"Tj) , we get the so-called 
Plemelj-Sokhotski formula 



P{z) 



I - 
— lim Im y 



I 



z — Xi + ir] 



(35) 



Via a change of basis, the sum on the right-hand side is 
equal to the trace of the matrix \{z + ir/)I — B] _1 , and 
hence p(z) is the limit where z goes to the real line of 
— (l/n7r)Im Tr(zl — B) _1 . In other words, the spectral 
density depends on the trace of the resolvent, and its 
average over the ensemble of model networks is given by 
the average of this quantity: 



1 



p(z) = ImTr^zI-B)- 1 ). 



(36) 



The normalized trace Tr(zl — B) _1 /n is called the Stielt- 
jes transform of B. 

The two most common ways to calculate the Stieltjes 
transform are either to expand the matrix (zl — B) _1 in 
powers of B and take the trace term by term, or to write 
the trace in terms of derivatives of a Fresnel integral and 
then employ the replica trick [40j . Here, however, we 
take a different approach inspired by work of Bai and 
Silverstein [l6|, |4l[ that allows us to calculate the average 
of the full resolvent. 

The resolvent is the inverse of a matrix whose off- 
diagonal elements are zero-mean random variables. Con- 
sider a general such matrix X and let us write it in terms 
of its first n — 1 rows and columns, plus the last row and 
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column, thus: 



X = 



X ra 



(37) 



Thus X„ is the matrix X with the nth row and column 
removed, and a is the nth column minus its last ele- 
ment x nn . 

Now consider the vector v = X _1 u, where u = 
(0, ...,0, 1). Let us break v into its first n — 1 ele- 
ments and its last element v = (vi|u n ), where clearly 
V n = [X _1 l . Then we have Xv = u and hence 

" L Inn 



X„Vi + v n a = 0, 



T 

a vi 



iV n = I- 



The first equation tells us that 

vi = -« n X~ 1 a, 
and substituting this result into the second gives 



TX" 1 ! = V n 

L inn " 



1 



a^X n ^a 



(38) 



(39) 



(40) 



To make further progress we assume that v n is nar- 
rowly peaked about its average value in the limit of large 
system size, meaning its variance about that value van- 
ishes as n becomes large. We will for the moment take 
this assumption as given, but it can be justified using 
results for concentration of measure of random quadratic 
forms , which apply provided vertex degrees are large 
(so that our results, like those of Section ing will be exact 
only for large degrees). 

If v n is narrowly peaked then the average of the re- 
ciprocal on the right-hand side of (|40|) is equal to the 
reciprocal of the average and 



- (a^X^a) 



(41) 



Furthermore, if v n is narrowly peaked then the average of 
Eq. (|39|) is (vi) = — u„(X„)(a) = since a is independent 
of X„ and (a) = 0. But the elements of Vi are equal 
to [X -1 ].; n and hence 



(42) 



for i =/= n. By the same method we can derive expressions 
for the inverse of X with any row and column i removed 
and hence show that 



i 



and 



<[ X-1 U=0 for^i- 



(43) 



(44) 



In other words, (X -1 ) is a diagonal matrix when n is 
large, with diagonal elements given by Eq. (|4*3"]) . 

But if this is true of X -1 , then by the same argument 
it must also be true of X" 1 . Hence, noting that a is 
independent of X.; , we have 



(a T X: 



jk 



(45) 

Returning now to Eq. (|36j) . the role of the matrix X 
in our problem is played by zl — B. As we noted ear- 
lier, the elements of the modularity matrix B (and hence 
also the elements of the vector a) have mean zero and 
variance k i kj/2m. Hence (a|) = fcjfcj/2m in Eq. (|4"5")) 
and 



(a r X- 1 a)=^([(zI-B. i) - 1 ] j3 ) 



1 3 



ii' 2m 



^Trp^I-B,)" 1 }], (46 ) 

where D is the diagonal matrix with elements ki and D; 
is the same matrix with the zth row and column removed. 

However, if Tr[D;((zI - B, i )~ 1 )]/2m tends to a well- 
defined limit as the network becomes large, then in this 
limit it must equal Tr[D((zI— B) -1 )] /2m — the omission, 
or not, of the ith row and column makes a vanishing 
difference for large n. Hence (|43| becomes 



<[(1-B)- 1 U = 



-fc 4 Tr[D((zI-B)-i)]/27 



(47) 



where we have made use of the fact that (Bn) = 0. At the 
same time, the off-diagonal elements of ({zl — B)" 1 ) arc 
zero by Eq. (|4"4")l , so the average of the resolvent matrix is 
diagonal, a result that will be crucial for several following 
developments. 

Without loss of generality, we now label the vertices of 
our network in order of increasing expected degree, and 
for convenience we define functions 7 z (a;) and k(x) of the 
continuous variable x thus: 

7 ,( i /n) = ([(zI-B)- 1 ] ij ), k(i/n) = ki. (48) 

Then for large n Eq. (|47[) becomes 

1 



[k{x)/c]$ Q k{y)~f z {y) Ay 



(49) 



where c = 2m /n is the average degree, as previously. 
The spectral density, Eq. (|36|) , is related to j z (x) by 



Pi*) 



lmg(z), 



where 



g(z)^-Tr((zI-B)- 1 )= f lz (x) dx, 
n Jo 



(50) 



(51) 
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which is just the ensemble average of the Stieltjes trans- 
form. To calculate g(z), we define the additional quantity 



Hz) 



-TrfD^zI-B)- 1 )] = - f k(x)j z (x)dx 
m cj 

k{x) dx 



2m 
c 



z — k(x)h(z) ' 



(52) 



where we have used Eq. (|49|) . Since we have labeled our 
vertices in order of increasing degree, k(x) is by definition 
the (nx)th-lowest degree in the network, or equivalently 
it is the functional inverse of the cumulative distribution 
function P(k) defined by 



P(k) 



p(k) dk' , 



(53) 



where p(k) is the expected degree distribution. Thus, 
changing variables from x to k, Eq. (|52j) can be written 



Hz) = - 

c 



fcdP(fc) 



Iq z — kh(z) ' 
or as either of the equivalent forms 

kp{k)dk f°° q(k)dk 



1 r 

h(z) = - \ 



kh{z) 



kh{z) 



(54) 



(55) 



where the (correctly normalized) probability distribution 



q(k) 



kp(k) 



(56) 



is known as the excess degree distribution in the net- 
works literature. This distribution, which arises often 
in the theory of networks, is the probability that the 
network vertex reached by following an edge has an ex- 
pected number k of edges attached to it other than the 
one we followed to reach the vertex. (The distribution 
looks slightly different from the form usually given [24| 
because it is expressed in terms of expected degree rather 
than actual degree.) 

If we can solve Eq. (|55|) for h(z) then we can calculate 
g(z) by substituting Eq. (|49|) into Eq. (f5"Tj) and again 
changing variables from x to k, to get 



9(z) 



p(k) dk 
z — kh(z) 



(57) 



This equation is similar in form to Eq. (|55[) , but note that 
it is the ordinary degree distribution p(k) that appears 
in the numerator, not the excess degree distribution. 

Alternatively, and more directly, we can calculate g(z) 
by multiplying both sides of (|4T))) by the right-hand de- 
nominator, integrating, and rearranging, to get 



l(z) = 



l + ch 2 {z) 



(58) 



Combining this result with Eq. (|50p now gives us the 
spectral density: 



lmh 2 {z), 



(59) 



where the imaginary part is, if necessary, calculated as 
the limit where z tends to the real line from above. 

Equations (|55| and ([59]) are precisely the equations, 
(fT9"| and ([23)1 . that we derived previously using the free 
convolution. 



V. SPECTRUM OF THE ADJACENCY MATRIX 

In the previous sections we have derived the spectral 
density of the modularity matrix. To calculate the cor- 
responding quantity for the adjacency matrix we make 
use of an argument of [IH, 0] as follows. The adjacency 
matrix can be written in terms of the modularity matrix 
as A = B + kk T /2m, where k is the n-element vector 
with elements fcj. Hence any eigenvalue/ vector pair z, v 
of the adjacency matrix satisfies 



B — ^ 

2m 



which can be rearranged to read 
k T v 



zv, 



2m 



(zI-B)" i k = v. 



(60) 



(61) 



Multiplying by k T , we then find that 

— k T (^I-B)- 1 k=l. 

2m 

Expanding k as a linear combination of the eigenvec- 
tors b, of B, this result can be written 



(62) 



V 



(k T b,) 2 



2m ' z - 

i—l 



fa 



(63) 



where fa are the eigenvalues of the modularity matrix. 

The solutions of this equation can be visualized as in 
Fig. [2] The solid curves represent the left-hand side of 
the equation, which has poles as shown at z = fa for 
all i. The dashed horizontal line represents the 1 on 
the right-hand side and the points at which it intercepts 
the curves are the solutions for z of (f6"3")> , which are the 
eigenvalues Ai of the adjacency matrix. If we number 
the eigenvalues of both A and B in order from largest to 
smallest, the geometry of Fig. [2] implies that the eigen- 
values must satisfy an interleaving condition of the form 
Ai > fa > A 2 > fa > ■ ■ ■ > K > Pn- In the limit of 
large n, where the spectral density of the modularity ma- 
trix becomes a smooth function and the eigenvalues are 
arbitrarily closely spaced, this implies that — > fa, so 
that asymptotically the spectral density of the adjacency 
matrix is the same as that of the modularity matrix. 

The only exception is the highest eigenvalue of the ad- 
jacency matrix Ai, which is bounded below by fa, but un- 
bounded above. To calculate this value we average (|6"2")l 
over the ensemble and recall, as demonstrated in Sec- 
tion IIV1 that ((zl — B) -1 ) is diagonal, and hence 



-Lk^((zI-B)- 1 )k=^ 

zm 



2m ^— u 



k*([(zI-B)-%). (64) 



9 





I 

n 


\ 




L 












\, \ , 




— 1 — h 


k 
\ 




\ / 


% y 


\ y 


i 



For the case of general degree distribution, we can 
use (|67p to eliminate h(z) in Eq. (f5"5j) to get 



1 



q(k) dk 



l-k/(z 2 -z)' 



(69) 



An exact solution to this equation requires us to perform 
the integral, but one can derive an approximate solution 
by expanding the denominator of the integrand: 



E 



q{k) dk, (70) 



FIG. 2: The solutions Xi to Eq. (|63[) correspond to the points 
where the left-hand side of the equation (solid curves) equals 1 
(dashed horizontal line). This implies that the values of the Xi 
are interleaved between the eigenvalues /3i of the modularity 
matrix. 



Combining this result with (p32"j) and using Eq. (|48p we 
then have 



1 



k 2 (x)^ z (x) dx = 1. 



(65) 



Taking Eq. (|49|) . multiplying by the right-hand denomi- 
nator and a further factor of k{x), then integrating over x, 
we get 



l r i 

2/ 



czh(z) — h(z) / k (x)7 z (x) dx = I k(x) dx. (66) 
Jo Jo 

And combining this result with Eq. (|65[) and noting that 
Jp k(x) dx = Jp 00 kp(k) dk = c, we have 



(z - l)h(z) = 1. 



(67) 



The solution of this equation for z gives us the leading 
eigenvalue Ai of the adjacency matrix. 

For the Poisson random graph, for example, this result, 
in combination with Eq. ([24]) . tells us that the leading 
eigenvalue takes the value c+1. This is not a new result — 
it is well known in the literature — but it is comforting to 
see that the formalism works. 

For the two-degree model of Eq. ([27} , we can use (|6"7| 
to eliminate h(z) from (f27|) and get 



Pidi + p 2 d 2 
(z-lf 



Pidi 



P2d 2 



z(z — 1) — d\ z(z — 1) — d-2 



(68) 



which gives us a cubic equation for z. For the parameter 
values used in Fig. Q] for example, d\ = 50, d 2 = 100, 
and pi = 1 — p 2 = j, we find that the leading eigenvalue 
of the adjacency matrix is z = 93.893 ... A numerical 
calculation for the same parameters is in good agreement, 
giving z = 93.896 ±0.017 for an average over 100 systems 
of size n = 10 000. 



1 00 



^ (z 2 - zY 

r—1 K y 



(71) 



where (. . .) q denotes an average over the excess degree 
distribution of Eq. (|56[) . If z 2 — z 3> fcmax, where fc max is 
the largest degree in the network, and noting that (k r ) q < 



fc ma XI w c have 



or 



+ 0[fc max /(z 2 -z)] 2 , 



(k 2 



(k) 



(72) 



(73) 



to leading order, where we have made use of (k) g = 
(k 2 )/(k). This result was derived previously by other 
means by Chung et al. [28|. 

Taking the example of the two degree model above 
again, this approximation gives 



Pid\+p 2 dl 
p\d\ +p 2 d 2 



(74) 



and for the parameter values of Fig. [T] we find that z ~ 
92.86, which differs by about 1% from the true value of 
93.89 given by Eq. 



VI. NETWORK HUBS 

The picture developed in the previous sections is one in 
which the spectrum of the adjacency matrix has two pri- 
mary components: a single leading eigenvalue plus a con- 
tinuous band of lower eigenvalues, which it shares with 
the modularity matrix. 

Let us examine more closely the continuous band and 
concentrate on the case of the modularity matrix, which 
is simpler since it has only the band and no separate lead- 
ing eigenvalue. Consider the eigenvalues that lie at the 
topmost edge of the band, which are the highest eigenval- 
ues of the modularity matrix. These eigenvalues are nor- 
mally associated with good bisections of the network into 
"communities" — if a good bisection exists then there will 
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be a corresponding high-lying eigenvalue whose eigenvec- 
tor's elements describe the split @. 

As we now argue, however, there is another mechanism 
that generates high-lying eigenvalues, namely the pres- 
ence of hubs in the network — vertices of unusually high 
degree — and the highest eigenvalues in the spectrum of 
the modularity matrix, and also the lowest, are often due 
to these hubs, while those corresponding to communities 
are somewhat smaller. As we will see, for hubs of suffi- 
ciently high degree, these eigenvalues can split off from 
the continuous band in a manner reminiscent of impu- 
rity states in condensed matter physics. In effect, the 
hub acts as an impurity in the network. 

To see how the addition of a hub to a network produces 
a high-lying eigenvalue, let the hub be vertex n and let 
B„ once again be the modularity matrix without the nth 
vertex (i.e., with the nth row and column removed), so 
that the full modularity matrix looks like this: 



B 



B, 



(75) 



Now, in an argument analogous to that of the pre- 
vious section, consider an eigenvector of this matrix 
v = (-vi\v n ). Then the eigenvector equation Bv = zv 
can be multiplied out to give the equations 

B„Vi + v n a = zvi, (76) 

a T vi + b nn v n = zv n . (77) 
The first of these can be rearranged to give 

v 1 = v n (zI-B n )- 1 a. (78) 

Then multiplying by a T and using the second equation 
gives 

a T (zI-B n )- 1 a = z-6„„. (79) 

Now we note that the ith element of a is an independent 
random variable with variance k n ki/2m and we can aver- 
age over the ensemble and apply Eq. (|45|) to rewrite the 
left-hand side, giving 



^Tr[D n <(,I-B„)-i>] 



(80) 



where D is the diagonal matrix with elements fc, as be- 
fore, D n is the same matrix with the nth row and col- 
umn removed, and we have made use of the fact that 
(bnn) = 0. We note, as previously, that if the quan- 
tity Tr(D„((zI — B„) _1 ))/2m tends to a limit as the 
network becomes large, then that limit is equal to the 
function h(z) defined in Eq. (f5"2"j) . Thus the eigenvalue z 
satisfies 
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FIG. 3: Graphical solution of Eq. (|81[) . The solid curves rep- 
resent the value of h(z) as a function of z, above and below the 
spectral band, and the hub eigenvalues, which are solutions 
of Eq. I)8ip . fall at the points where this curve intersects the 
straight line z/k n , represented by the dashed diagonal line. 
The slope of this line is l/k n , and hence when k n is large 
enough the lines intersect — case (a) — and we have two hub 
eigenvalues, one above and one below the band (marked by 
dots). Case (b) is the borderline case. If k n is any less than 
this value then there is no intersection and the highest and 
lowest eigenvalues will be those at the band edges. 

Substituting this expression into Eq. (|55|) and rearrang- 
ing, we get an explicit expression for the eigenvalue thus: 



kr. 



kp(k) dk 



(82) 



This calculation also extends to the case where there is 
more than one hub in the network. Because the hub 
is treated no differently from any other network vertex, 
the same arguments apply if we add a second hub, or 
more, after the first. Equation (|82|) will give the correct 
eigenvalue for each hub separately. 

Once again, our ability to actually solve for the value 
of z will depend on whether we can do the integral in 
Eq. (|82j) (although one could also evaluate the integral 
numerically) . In the special case where the hub degree k n 
is much larger than the expected degree of any of the 
other vertices, so that k n — k ~ k n in the denominator of 
the integrand, the expression simplifies to 



c 



poo 

/ kp(k) dk = k r - 
Jo 



(83) 



Kz) 



(81) 



and hence z = \fk~^. 

The solutions of Eq. (fSTj) can be represented graphi- 
cally as in Fig. [3] The curves in the figure represent the 
function h(z) and the diagonal lines represent z/k n . The 
point where the two cross give the eigenvalues. As the fig- 
ure shows, when the expected degree k n of the nth vertex 
is large enough, the equation has two solutions, one for 
low z and one for high and both given by Eq. (|52")h that 
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are separate from the continuous spectrum of eigenvalues 
we calculated in Section HIT1 

How high a degree does a hub have to have to gen- 
erate eigenvalues of this kind? The answer can be seen 
from Fig. [3] — k n must be large enough for the line z/k n 
to intercept the curve of h(z). Thus there is a critical 
value of k n , represented by the steeper diagonal in the 
figure, below which the hub eigenvalues vanish. Below 
this point, the highest eigenvalue will fall at the edge 
of the continuous band as normal and there will be no 
special hub eigenvalue. We can derive an expression for 
the transition point by observing that, as shown in Sec- 
tion [BIB] the slope of h(z) diverges at the band edge, 
which implies that dz/dk n = 0. Differentiating Eq. (|82|) 
and setting the result to zero, we find that the critical 
value of k n is the solution of 



30 - 



kp{k) dk 
kn k 



k 2 p(k) dk 

{k n 



(84) 



For example, in the case of the Poisson random graph 
this implies that the transition takes place at the point 
where c/(k n — c) = c 2 /(k n — c) 2 , i.e., when k n = 2c. Thus 
we must have k„ > 2c for the hub to have an effect on 
the spectrum. 

This gives us a working definition of what we mean by 
a "hub" in a network. It depends, not surprisingly, on 
the degree distribution of the rest of the network — what 
it takes to stand out in a crowd depends on the rest of the 
crowd. But in the Poisson random graph, for instance, 
a hub is a hub, in spectral terms, if its degree is greater 
than twice the average in the rest of the network. This is 
a somewhat surprising result, given that vertices of high 
degree are easily spotted long before this point is reached, 
at least for large c. Since the standard deviation of the 
degree distribution is -y/c, a vertex with degree twice the 
mean is yfc standard deviations above the mean, which 
is a large number for large c. 

Nonetheless, the result does appear to be correct. Fig- 
ure |4] shows the results of numerical calculations of the 
largest eigenvalue of the modularity matrix for a Poisson 
random graph with a single additional hub of expected 
degree k n , as a function of k n . As the figure shows, the 
eigenvalue obeys Eq. ([82]) quite closely until k n falls be- 
low 2c (the vertical dashed line). Past this point, the 
leading eigenvalue assumes the same value 2-y/c as in a 
standard Poisson random graph with no hub (the hori- 
zontal line), even though the hub may still be present. 

Putting together our principal observations, we have 
now developed quite a complete picture of the spectrum 
of the configuration model. We expect the spectrum to 
have two main parts, plus a third when the degree dis- 
tribution implies the presence of hubs: 

1 . There is a single eigenvalue given by the solution of 
Eq. (p7|) , which will normally be the leading eigen- 
value. 
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FIG. 4: The largest eigenvalue of the modularity matrix for 
a Poisson random graph of mean degree c = 100 plus a single 
additional hub of expected degree k n . Points are numerical 
results, averaged over 1000 networks of n = 10 000 vertices 
each. Statistical errors on the measurements are smaller than 
the points in all cases. The solid curve is Eq. (|82[) . which gives 
z — kn/Vkn — c in this case, and the horizontal dashed line 
represents the value z = 2 v / c = 20, which is the lower limit 
on the eigenvalue set by the edge of the continuous spectral 
band. The vertical dashed line represents the critical value 
k n = 200 of the hub degree, set by Eq. ((84)) . 



bounded, both above and below, and have edges 
that decay to zero as a square root. 

3. If there are hubs in the network, then there will 
be additional eigenvalues outside the band at both 
ends, given by Eq. (|82[). Each hub contributes two 
eigenvalues, one at each end of the band. 



A. Localization around hubs 

One can also look at the eigenvector corresponding to 
a hub eigenvalue, which turns out to be heavily localized 
around the hub vertex. All the elements of the eigenvec- 
tor, except for the element v n corresponding to the hub 
itself, are given in terms of v n by Eq. (|78[) . For given a, 
the expected value of the ith component is 



((zI-B„)- 1 )a 



((^I-B,,)- 1 ) 



(85) 



2. There is a continuous band, given by Eq. ([23]) . For 
bounded degree distributions the band will also be 



where we have once again made use of the fact that ((zl— 
Bn)- 1 ) is diagonal (see Eq. (|4l)0 . 

The ith element of the vector a takes the value a, = 
1 — kik n /2m for vertices i that are connected to the hub 
and —kik n /2m for those that are not. Hence, in the 
limit of large n, eigenvector elements corresponding to 
neighbors of the hub will be of order a constant, with 
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expected value 

Vi = v n j z (i/n) 



z — kih(z) z(l — ki/k n ) 



(86) 



with z given by Eq. (|82[) . while the remaining elements 
will be of order 1/n. 

The value of v n can be determined by insisting that 
the complete eigenvector be normalized. Using Eq. ([75)1 
we can write the normalization condition in the form 



1 



^ + |v 1 | 2 =^[l + a r (.I-B)- 2 a] 



-ia^zI-B)-^ 

dz 



(87) 



When we average over the ensemble we have, by analogy 
with Eq. (|46|) . 



<a T (zI - B)-^) = ^ Tr(D(zI - B)" 1 ) = k n h(z), 



and hence ([87]) implies that 



1 



l-k n h'(z)' 



(88) 



(89) 



where h'(z) denotes the first derivative of h(z), and we 
are assuming once again that the vector element v n is 
narrowly peaked about its expected value. Note that 
h'(z) is negative at both the positive and negative band 
edges, and diverges to -co as we approach the band edge. 
Thus v n — > as we approach the transition at the which 
the hub eigenvalue disappears. 

The results above apply to the hub eigenvectors at both 
ends of the spectral band, there being two eigenvalues for 
each hub vertex, one at either end, as shown in the previ- 
ous section. Both eigenvectors will have a single element 
of order 1 in the position corresponding to the hub itself, 
elements of order 1/z in the positions corresponding the 
neighbors of the hub (see Eq. (j86f ). and all other elements 
of order 1/n. In other words, the both eigenvectors are 
strongly localized in the neighborhood of the hub. The 
only qualitative difference between the two eigenvectors 
is in the sign of the elements corresponding to the neigh- 
bors which, because of Eq. ([SB]) , will have the same sign 
as v n for the positive eigenvalue and the opposite sign for 
the negative one. 

As an example, consider again the Poisson random 
graph, for which h(z) is given by Eq. (f!H|) and z is 
given by Eq. (|82j) to be ±k n /y/k n — c, so that h'(z) = 
— 1/ (k n — 2c) and the expected values of the eigenvector 
elements at both ends of the spectrum satisfy 



2 k n 



c) / [k n — c) for i = n, 

c) / (k n — c) 2 for i a neighbor of n, 

otherwise, 

(90) 

in the limit of large network size. Figure [5] shows a com- 
parison of these predictions with numerical results for 
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FIG. 5: Values of elements of the leading eigenvector of 
the modularity matrix for a Poisson random graph with 
n — 10 000 vertices and mean degree c = 100, with a single 
added hub of degree k n . Main figure: value of the vector ele- 
ment corresponding to the hub itself. Inset: average value of 
the elements corresponding to the hub's immediate network 
neighbors. Points are numerical results, averaged over 100 
networks each; curves are the analytic prediction, Eq. (|90|) . 
Statistical errors are smaller than the data points in all cases. 



actual networks. As the figure shows, the agreement is 
once again good, although, as with some of the other 
calculations, there are small disparities close to the tran- 
sition at which the hub eigenvalue meets the band edge 
(which is at k n = 200 in this case). 



VII. CONCLUSIONS 

In this paper we have studied the spectra of the adja- 
cency and modularity matrices of random networks with 
given expected degrees. Our principal findings are that 
the spectral densities of the adjacency and modularity 
matrices are the same in the limit of large system size, 
except that the adjacency matrix has an additional high- 
est eigenvalue, and that the spectral densities are given 
by the free multiplicative convolution of the degree dis- 
tribution with a Wigner semicircle distribution. We have 
confirmed these results with numerical studies of actual 
networks generated according to the model. The spectra 
show strong departures from the classical semicircle law, 
in agreement with numerical studies by previous authors. 

We have also studied the effect of network hubs, ver- 
tices of unusually high degree, and find that when their 
degree is sufficiently large these give rise to eigenvalues 
outside the main band of the spectrum, akin to impurity 
states in condensed matter systems. We have derived an 
explicit formula for these hub eigenvalues and we show 
that the corresponding eigenvectors are strongly localized 
around the hubs themselves. 

In addition to their relevance to partitioning, commu- 
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nity structure, and dynamical systems on networks, the 
techniques developed here could form a starting point for 
spectral calculations in more elaborate networks. There 
has, for instance, been recent interest in the spectral 
properties of community structured networks (45l l46j . 
but calculations have been limited to models with Pois- 
son degree distribution. Applications of the methods pre- 
sented here to such networks could lead to new results for 
structured networks with nontrivial degree distributions. 



Acknowledgments 

The authors thank Lenka Zdeborova for useful con- 
versations. This work was funded in part by the Na- 
tional Science Foundation under grants CCF-1116115 
and DMS-1 107796 and by the Army Research Office un- 
der MURI grant W911NF-11-1-0391. 



M. E. J. Newman, The structure and function of complex [20 

networks. SI AM Review 45, 167-256 (2003). 

S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.- 

U. Hwang, Complex networks: Structure and dynamics. [21 

Physics Reports 424, 175-308 (2006). 

M. Fiedler, Algebraic connectivity of graphs. Czech. 

Math. J. 23, 298-305 (1973). [22 

A. Pothen, H. Simon, and K.-P. Liou, Partitioning sparse 
matrices with eigenvectors of graphs. SIAM J. Matrix 
Anal. Appl. 11, 430-452 (1990). [23 

B. Bollobas, C. Borgs, J. Chayes, and O. Riordan, Per- 
colation on dense graph sequences. Annals of Probability 
38, 150-183 (2010). [24 
M. E. J. Newman, Modularity and community structure 

in networks. Proc. Natl. Acad. Set. USA 103, 8577-8582 
(2006). [25 
S. Fortunato, Community detection in graphs. Phys. Rep. 
486, 75-174 (2010). 
A. Arenas, A. Dfaz-Guilera, J. Kurths, Y. Moreno, and [26 

C. Zhou, Synchronization in complex networks. Physics 
Reports 469, 93-153 (2008). 

A. Barrat, M. Barthelemy, and A. Vespignani, Dynamical 
Processes on Complex Networks. Cambridge University [27 
Press, Cambridge (2008). 

I. J. Farkas, I. Derenyi, A.-L. Barabasi, and T. Vicsek, 

Spectra of "real-world" graphs: Beyond the semicircle [28 

law. Phys. Rev. E 64, 026704 (2001). 

R. Solomonoff and A. Rapoport, Connectivity of random 

nets. Bulletin of Mathematical Biophysics 13, 107-117 [29 

(1951). 

P. Erdos and A. Renyi, On the evolution of random 
graphs. Publications of the Mathematical Institute of the [30 
Hungarian Academy of Sciences 5, 17-61 (1960). 
L. Arnold, On Wigner's semicircle law for the eigenval- 
ues of random matrices. Probability Theory and Related [31 
Fields 19, 191-198 (1971). 

Z. Fiiredi and J. Komlos, The eigenvalues of random sym- [32 
metric matrices. Combmatorica 1, 233-241 (1981). 
G. Anderson, A. Guionnet, and O. Zeitouni, An Introduc- 
tion to Random Matrices. Cambridge University Press, [33 
Cambridge (2010). 

Z. Bai and J. Silverstein, Spectral Analysis of Large Di- 
mensional Random Matrices. Springer, Berlin (2010). [34 
P. van Mieghem, Graph Spectra for Complex Networks. 
Cambridge University Press, Cambridge (2011). 
L. Erdos, A. Knowles, H. Yau, and J. Yin, Spectral 
statistics of Erdos-Renyi graphs I: Local semicircle law. [35 
Preprint \arXiv:1103.1919\ 

T. Tao, Topics in Random Matrix Theory. American 
Mathematical Society, Providence, RI (2012). [36 



T. Tao and V. Vu, Random matrices: The uni- 
versality phenomenon for Wigner ensembles. Preprint 
\arXiv:1202.0068\ 

B. Bollobas, A probabilistic proof of an asymptotic for- 
mula for the number of labelled regular graphs. European 
Journal of Combinatorics 1, 311-316 (1980). 
M. Molloy and B. Reed, A critical point for random 
graphs with a given degree sequence. Random Structures 
and Algorithms 6, 161-179 (1995). 

M. Molloy and B. Reed, The size of the giant component 
of a random graph with a given degree sequence. Combi- 
natorics, Probability and Computing 7, 295-305 (1998). 
M. E. J. Newman, S. H. Strogatz, and D. J. Watts, Ran- 
dom graphs with arbitrary degree distributions and their 
applications. Phys. Rev. E 64, 026118 (2001). 
R. Cohen, K. Erez, D. ben-Avraham, and S. Havlin, Re- 
silience of the Internet to random breakdowns. Phys. Rev. 
Lett. 85, 4626-4628 (2000). 

D. S. Callaway, M. E. J. Newman, S. H. Strogatz, and 
D. J. Watts, Network robustness and fragility: Percola- 
tion on random graphs. Phys. Rev. Lett. 85, 5468-5471 
(2000). 

S. N. Dorogovtsev, A. V. Goltsev, J. F. F. Mendes, and 
A. N. Samukhin, Spectra of complex networks. Phys. 
Rev. E 68, 046109 (2003). 

F. Chung, L. Lu, and V. Vu, Spectra of random graphs 
with given expected degrees. Proc. Natl. Acad. Sci. USA 
100, 6313-6318 (2003). 

F. Chung and L. Lu, The average distances in random 
graphs with given expected degrees. Proc. Natl. Acad. 
Sci. USA 99, 15879-15882 (2002). 

D. Voiculescu, K. Dykema, and A. Nica, Free Random 
Variables, Volume 1. American Mathematical Society, 
Providence, RI (1992). 

G. J. Rodgers and A. J. Bray, Density of states of a sparse 
random matrix. Phys. Rev. B 37, 3557-3562 (1988). 

N. Rao and A. Edelman, The polynomial method for 
random matrices. Foundations of Computational Mathe- 
matics 8, 649-702 (2008). 

S. Olver and R. Nadakuditi, Numerical computation 
of convolutions in free probability theory. Preprint 
\arXiv:1203.1958\ 

N. Rao and R. Speicher, Multiplication of free ran- 
dom variables and the S-transform: The case of vanish- 
ing mean. Electronic Communications in Probability 12, 
248-258 (2007). 

S. Molchanov, L. Pastur, and A. Khorunzhii, Limiting 
eigenvalue distribution for band random matrices. Theo- 
retical and Mathematical Physics 90, 108-118 (1992). 
D. Shlyakhtenko, Random Gaussian band matrices and 



14 



freeness with amalgamation. International Mathematics 
Research Notices 20, 1013-1025 (1996). 
[37] G. Anderson and O. Zeitouni, A CLT for a band matrix 
model. Probability Theory and Related Fields 134, 283- 
338 (2006). 

[38] G. Casati and V. Girko, Wigner's semicircle law for 
band random matrices. Random Operators and Stochas- 
tic Equations 1, 15-22 (2009). 

[39] Z. Bai and L. Zhang, The limiting spectral distribution 
of the product of the Wigner matrix and a nonnegative 
definite matrix. Journal of Multivariate Analysis 101, 
1927-1949 (2010). 

[40] S. F. Edwards and R. C. Jones, The eigenvalue spectrum 
of a large symmetric random matrix. J. Phys. A 9, 1595- 
1603 (1976). 

[41] Z. Bai, Methodologies in spectral analysis of large- 
dimensional random matrices: A review. Statist. Smica 
9, 611-677 (1999). 

[42] N. Karoui and H. Koesters, Geometric sensitivity of 



random matrix results: Consequences for shrinkage es- 
timators of covariance and related statistical methods. 
Preprint arXiv: 1 1 05. 1404] 

[43] M. Capitaine, C. Donati-Martin, and D. Feral, The 
largest eigenvalues of finite rank deformation of large 
Wigner matrices: Convergence and nonuniversality of the 
fluctuations. Annals of Probability 37, 1-47 (2009). 

[44] F. Benaych-Georges and R. R. Nadakuditi, The eigenval- 
ues and eigenvectors of finite, low rank perturbations of 
large random matrices. Advances in Mathematics 227, 
494-521 (2011). 

[45] S. Chauhan, M. Girvan, and E. Ott, Spectral properties 
of networks with community structure. Phys. Rev. E 80, 
056114 (2009). 

[46] R. R. Nadakuditi and M. E. J. Newman, Graph spectra 
and the detectability of community structure in networks. 
Phys. Rev. Lett. 108, 188701 (2012). 



